# Section 1 工作区设置和R包加载 ----------------------------------------------------

setwd("E://OneDrive//研究//8-学位论文//数据分析处理//第二版")

library(xlsx) # Read, Write, Format Excel 2007 and Excel 97/2000/XP/2003 Files
library(plyr) # Tools for Splitting, Applying and Combining Data
library(tidyr) # Tidy Messy Data
library(agricolae) # Statistical Procedures for Agricultural Research
library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics
library(gridExtra) # Miscellaneous Functions for "Grid" Graphics
library(ggpubr) # 'ggplot2' Based Publication Ready Plots


# Section 2 自定义参数 ----------------------------------------------------
rm(list=ls())

windowsFonts(myFont = windowsFont("微软雅黑")) #绘图字体设置，会报Warning，但不影响出图
color_yb2 <- c("#1f78b4", "#ff7f00") #双色色号
color_yb4 <- c("#1f78b4", "#ff7f00", "#7ab5e2", "#ffd92f") #四色色号

# Section 3 自定义函数 ---------------------------------------------------------

bootstrap.yb <- function(data, tv, cv, group_c = "Target_Type", group_t, n_iter, n_samp) {  
  library(plyr) # Tools for Splitting, Applying and Combining Data                                #dlply拆分数据框
  library(data.table)                          #data.table包tstrsplit、rbindlist函数
  
  group <- c(group_t, group_c)                 #分组变量
  
  data_c <- dlply(data, group_c)               #按组将数据拆分进list
  C_names <- names(data_c)                     #将分组组合
  C_g <- as.data.frame(tstrsplit(C_names, ".", fixed=TRUE))
  #分组组合拆为数据框
  names(C_g) <- group_c                        #命名分组数据框列名
  
  
  data_t <- dlply(data, group)                 #按组将数据拆分进list
  T_names <- names(data_t)                     #将分组组合
  T_g <- as.data.frame(tstrsplit(T_names, ".", fixed=TRUE))
  #分组组合拆为数据框
  names(T_g) <- group                          #命名分组数据框列名
  
  
  E <- list()                                  #效应list
  for (i in c(1:length(data_t))) {             #按分组数循环
    bs <- c()
    for (j in c(1:n_iter)) {
      set.seed(j)
      bs <- append(bs, mean(sample(data_t[[i]][tv][,1], n_samp, replace = TRUE), na.rm = T) - mean(sample(data_t[[i]][cv][,1], n_samp, replace = TRUE), na.rm = T))
    }
    E[[i]] <- cbind(data.frame(T_g[i,]), bs, row.names = NULL)
    names(E[[i]]) <- c(group, "E")             #给数据框命名
  }
  return(E)
}
#data：需数据框格式, tv：处理变量, cv：对照变量, group_c：对照分组（默认物种）, group_t：处理分组, n_iter：迭代次数, n_samp：抽样数

mq.yb <- function(data, v, group) {               
  library(plyr) # Tools for Splitting, Applying and Combining Data                                #dlply拆分数据框
  datac <- ddply(data, group,
                 .fun = function(xx, col) {
                   c(mean = mean     (xx[[col]], na.rm=T),
                     down = quantile (xx[[col]], 0.025, na.rm=T),
                     up   = quantile (xx[[col]], 0.975, na.rm=T),
                     p    = length(which(xx[[col]]>0))/length(xx[[col]])
                   )
                 },
                 v
  )
  names(datac) <- c(group, "mean", "down", "up", "p")
  datac <- datac[which(is.nan(datac$mean) == F),]
  return(datac)
}#将data数据框按group分组计算v的平均值mean，2.5%和97.5%的分位数down、up，和大于0的比例p

replace.yb <- function(vector, from, to) {       
  x <- vector
  for (i in c(1:length(from))) {
    x[x == from[i]] <- to[i]
  }
  return(x)
}#批量替换向量vector中的元素，将from元素替换为to元素

inter.yb <- function(list, ab, a, b, v) {      
  out <- list()
  for (i in c(1:length(ab))) {
    up <- max(quantile(list[[a[i]]][v], 0.975, na.rm=T), quantile(list[[b[i]]][v], 0.975, na.rm=T))
    down <- min(quantile(list[[a[i]]][v], 0.025, na.rm=T), quantile(list[[b[i]]][v], 0.025, na.rm=T))
    E <- as.vector(unlist(list[[a[i]]][v] + list[[b[i]]][v]))
    Emean <- mean(E, na.rm=T)
    Edown <- quantile (E, 0.025, na.rm=T)
    Eup   <- quantile (E, 0.975, na.rm=T)
    list[[ab[i]]]$inter[unlist(list[[ab[i]]][v]) - E > Edown - Emean & unlist(list[[ab[i]]][v]) - E < Eup - Emean] <- "Additive"
    if (mean(unlist(list[[a[i]]][v]), na.rm=T) / mean(unlist(list[[b[i]]][v]), na.rm=T) < 0) {
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean > 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) > 0] <- "+Synergistic"
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean > 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) < 0] <- "-Antagonistic"
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean < 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) < 0] <- "+Antagonistic"
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean < 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) > 0] <- "-Synergistic"
    } else {
      list[[ab[i]]]$inter[list[[ab[i]]][v] - E < Edown - Emean] <- "Antagonistic"
      list[[ab[i]]]$inter[list[[ab[i]]][v] - E > Eup - Emean] <- "Synergistic"
    }
    out[[i]] <- list[[ab[i]]]
  }
  return(out)
}#计算交互效应，list列表，ab、a、b交互、因子1、因子2序号向量，v比较的数据

relative.yb <- function(data, Col_DV, Col_IV, Col_group = NULL) {
  df <- data
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  data_list <- list()
  out_list <- list()
  out <- data.frame()
  for (i in c(1:length(group))) {
    data_list[[i]] <- df[df$group == group[i],]
    out_list[[i]] <- df[df$group == group[i], c(Col_IV, Col_group)]
    for (j in c(1:length(Col_IV))) {
      if(j == 1) {
        CK_row <- which(data_list[[i]][,Col_IV[j]] == 0)
      } else {
        CK_row <- intersect(CK_row, which(data_list[[i]][,Col_IV[j]] == 0))
      }
    }
    for (j in c(1:length(Col_DV))) {
      x0 <- mean(data_list[[i]][CK_row,Col_DV[j]])
      relative <- (data_list[[i]][,Col_DV[j]] - x0)/x0 *100
      out_list[[i]] <- cbind(out_list[[i]], relative)
      names(out_list[[i]])[names(out_list[[i]]) == "relative"] <- Col_DV[j]
    }
    out <- rbind(out, out_list[[i]])
  }
  return(out)
}

ANOVA.yb <- function(df, Col_DV, Col_IV, Col_group = NULL) {
  library(tidyr) # Tidy Messy Data # Tidy Messy Data # Tidy Messy Data
  library(car) # Companion to Applied Regression
  library(plyr) # Tools for Splitting, Applying and Combining Data
  library(dplyr) # A Grammar of Data Manipulation
  library(agricolae) # Statistical Procedures for Agricultural Research
  
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  
  result <- list()
  out <- list()
  
  for (i in c(1:length(group))) {
    result[[i]] <- list()
    out[[i]] <- list()
    df1 <- df[df$group == group[i],]
    df1 <- unite(df1, "IV", all_of(Col_IV), sep = "_", remove = F)
    df1$IV <- paste("IV", df1$IV, sep = "_")
    df1$IV <- as.factor(df1$IV)
    IV <- unique(df1$IV)
    
    for (j in c(1:length(Col_DV))) {
      result[[i]][[j]] <- list()
      length2 <- function (x, na.rm=T) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
      }
      datac <- ddply(df1, c(Col_IV, "IV"), .drop=T,
                     .fun = function(xx, col) {
                       c(N    = length2 (xx[[col]], na.rm=T),
                         mean = mean    (xx[[col]], na.rm=T),
                         sd   = sd      (xx[[col]], na.rm=T),
                         min  = min     (xx[[col]], na.rm=T),
                         max  = max     (xx[[col]], na.rm=T)
                       )
                     },
                     Col_DV[j]
      )
      datac$se <- datac$sd / sqrt(datac$N) 
      datac$down <- datac$mean - datac$se
      datac$up <- datac$mean + datac$se
      result[[i]][[j]][[3]] <- datac
      
      model <- paste(Col_DV[j], paste(Col_IV, collapse = "*"), sep = "~")
      normal <- data.frame(group = "", DV = "", IV = "", KS = "", KS_P = "", SW = "", SW_P = "")[-1,]
      for (k in c(1:length(IV))) {
        sample <- df1[df1$IV == IV[k], Col_DV[j]]
        ks <- ks.test(sample, pnorm, mean(sample), sd(sample))
        if (sd(sample) != 0) {
          sw <- shapiro.test(sample)
        } else {
          sw <- list()
          sw[[1]] <- 0
          sw[[2]] <- 0
        }
        normal <- rbind(normal, c(as.character(group[i]), as.character(Col_DV[j]), as.character(IV[k]), 
                                  round(ks[[1]], 4), round(ks[[2]], 4), 
                                  round(sw[[1]], 4), round(sw[[2]], 4)))
      }
      names(normal) <- c("group", "DV", "IV", "KS", "KS_P", "SW", "SW_P")
      result[[i]][[j]][[1]] <- as.data.frame(leveneTest(eval(parse(text = model)), data = df1, center = mean))
      result[[i]][[j]][[1]] <- cbind(rep(group[i], length(result[[i]][[j]][[1]][,1])),
                                     rep(Col_DV[j], length(result[[i]][[j]][[1]][,1])),
                                     result[[i]][[j]][[1]])
      names(result[[i]][[j]][[1]])[1:2] <- c("group", "DV")
      result[[i]][[j]][[2]] <- as.data.frame(summary(aov(eval(parse(text = model)), data = df1))[[1]])
      result[[i]][[j]][[2]] <- cbind(rep(group[i], length(result[[i]][[j]][[2]][,1])),
                                     rep(Col_DV[j], length(result[[i]][[j]][[2]][,1])),
                                     row.names(result[[i]][[j]][[2]]),
                                     result[[i]][[j]][[2]])
      names(result[[i]][[j]][[2]])[1:3] <- c("group", "DV", "source")
      SNK <- SNK.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      SNK <- data.frame(IV = rownames(SNK), SNK = SNK$groups)
      LSD <- LSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV", p.adj="none")$group
      LSD <- data.frame(IV = rownames(LSD), LSD = LSD$groups)
      duncan <- duncan.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      duncan <- data.frame(IV = rownames(duncan), duncan = duncan$groups)
      HSD <- HSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      HSD <- data.frame(IV = rownames(HSD), HSD = HSD$groups)
      mc <- merge(merge(SNK, LSD, by = "IV"), merge(duncan, HSD, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- merge(result[[i]][[j]][[3]], merge(normal, mc, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- result[[i]][[j]][[3]][, c("group", "DV", Col_IV, "N", "mean", "sd", "min", "max", "se", "down", "up", "KS", "KS_P", "SW", "SW_P", "SNK", "LSD", "duncan" ,"HSD")]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$SNK) > 2,"SNK"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$SNK)[nchar(result[[i]][[j]][[3]]$SNK) > 2]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$LSD) > 2,"LSD"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$LSD)[nchar(result[[i]][[j]][[3]]$LSD) > 2]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$duncan) > 2,"duncan"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$duncan)[nchar(result[[i]][[j]][[3]]$duncan) > 2]
      result[[i]][[j]][[3]][nchar(result[[i]][[j]][[3]]$HSD) > 2,"HSD"] <- sub('^(.).*(.)$', '\\1-\\2', result[[i]][[j]][[3]]$HSD)[nchar(result[[i]][[j]][[3]]$HSD) > 2]
    }
    for (j in c(1:length(Col_DV))) {
      if (j == 1) {
        out[[i]][[1]] <- result[[i]][[j]][[1]]
        out[[i]][[2]] <- result[[i]][[j]][[2]]
        out[[i]][[3]] <- result[[i]][[j]][[3]]
      } else {
        out[[i]][[1]] <- rbind(out[[i]][[1]], result[[i]][[j]][[1]])
        out[[i]][[2]] <- rbind(out[[i]][[2]], result[[i]][[j]][[2]])
        out[[i]][[3]] <- rbind(out[[i]][[3]], result[[i]][[j]][[3]])
      }
    }
  }
  return(out)
} #方差分析；输入(列名)：df-数据表(数据框)、DV-因变量(向量)、VI-自变量/固定因子(向量)、group-多重比较分组(向量)
#输入：result[[group]][[DV]][[1-levene(n>50-KS); 2-anova; 3-summary & normality_test & multiple_comparisons]]

analysis.yb <- function(data, Col_DV, Col_IV, Col_group, Col_species, species) {
  library(tidyr) # Tidy Messy Data # Tidy Messy Data # Tidy Messy Data
  df <- data
  df <- unite(df, "species", all_of(Col_species), sep = "_", remove = F)
  df$species <- factor(df$species, levels = species)
  #species <- unique(df$species)
  df1 <- list()   #df1未经转换的数据（[[1]]为全部数据，[[2]]之后为按物种分类后的相对值）
  df1[[1]] <- df[,c(Col_species, Col_group, Col_IV, Col_DV)]
  for (i in c(1:length(species))) {
    df_species <- df[df[,Col_species] == species[i], c(Col_group, Col_IV, Col_DV)]
    df1[[i+1]] <- cbind(rep(species[i], length(df_species[,1])),
                        relative.yb(data = df_species, Col_DV, Col_IV, Col_group))
    names(df1[[i+1]])[1] <- Col_species
  }
  df2 <- list()
  for (i in c(1:length(df1))) {
    if (i == 1) {
      df2[[i]] <- ANOVA.yb(df = df1[[i]], Col_DV, Col_IV = c(Col_IV, Col_group, Col_species))
    } else {
      df2[[i]] <- ANOVA.yb(df = df1[[i]], Col_DV, Col_IV = c(Col_IV, Col_group))
    }
  }
  return(df2)
} #只用于博士毕业论文，未做兼容性设计 输出1原始值方差分析，2以后不同物种的方差分析

fig.yb <- function(data, Col_IV, x1, x2, Col_DV = "DV", y1, y2, 
                   Col_group = "group", g1, g2, pointrange = T) {
  library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics # Create Elegant Data Visualisations Using the Grammar of Graphics
  library(tidyr) # Tidy Messy Data # Tidy Messy Data # Tidy Messy Data
  
  data <- unite(data, "Group", all_of(Col_group), sep = "_", remove = F)
  for (i in c(1:length(g1))) {
    data$Group[data$Group == g1[i]] <- g2[i]
  }
  data$Group <- factor(data$Group, levels = g2)
  
  data <- unite(data, "IV", all_of(Col_IV), sep = "_", remove = F)
  data$IV <- paste("IV", data$IV, sep = "_")
  for (i in c(1:length(x1))) {
    data$IV[data$IV == x1[i]] <- x2[i]
  }
  data$IV <- factor(data$IV, levels = x2)
  
  Fig <- list()
  
  if (pointrange) {
    for (i in c(1:length(y1))) {
      Fig[[i]] <- ggplot(data = data[data[,Col_DV] == y1[i],])+
        geom_hline(yintercept = 0, linetype = "dashed", color = "#727272")+
        geom_pointrange(aes(x = IV, y = mean, ymin = down, ymax = up, color = Group),
                        position = position_dodge(.75), size = 1)+
        geom_text(aes(x = IV, y = up, label = duncan), color = "#000000", angle = 90,
                  position = position_dodge2(.75), hjust = -0.5, size = 4.5, fontface = "bold", show.legend = F)+
        scale_color_manual(values = color_yb2)+
        scale_fill_manual(values = color_yb2)+
        scale_y_continuous(labels=function(x) paste0(x,"%"),
                           limits = c(min(data[data[,Col_DV] == y1[i],]$down), 
                                      1.2*max(data[data[,Col_DV] == y1[i],]$up) - 0.2*min(data[data[,Col_DV] == y1[i],]$down)))+
        labs(y = y2[i])+
        theme_bw()+
        theme(legend.position="top",
              legend.title=element_blank(),
              legend.text = element_text(size = 12),
              axis.text.x = element_text(size = 12, face = "bold"),
              axis.title.x = element_blank(),
              axis.title.y = element_text(size = 14),
              axis.text.y = element_text(size = 12),
              text = element_text(family = "myFont"))
    }
  } else {
    for (i in c(1:length(y1))) {
      Fig[[i]] <- ggplot(data = data[data[,Col_DV] == y1[i],])+
        geom_col(aes(x = IV, y = mean, color = Group, fill = Group), width = 0.75, position = position_dodge())+
        geom_errorbar(aes(x = IV, ymin = down, ymax = up, color = Group), 
                      width = 0.2, position = position_dodge(.75), size = 1)+
        geom_text(aes(x = IV, y = up, label = duncan, group = Group), color = "#000000", angle = 90,
                  position = position_dodge(.75), hjust = -0.5, size = 4.5, fontface  = "bold", show.legend = F)+
        scale_color_manual(values = color_yb4)+
        scale_fill_manual(values = color_yb4)+
        scale_y_continuous(expand = c(0,0),
                           limits = c(0, max(data[data[,Col_DV] == y1[i],]$up)*1.2))+
        labs(y = y2[i])+
        theme_bw()+
        theme(legend.position="top",
              legend.title=element_blank(),
              legend.text = element_text(size = 12),
              axis.text.x = element_text(size = 12, face = "bold"),
              axis.title.x = element_blank(),
              axis.title.y = element_text(size = 14),
              axis.text.y = element_text(size = 12),
              text = element_text(family = "myFont"))
    }
  }
  return(Fig)
} #博士毕业论文绘图

temperature.yb <- function(data) {
  library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics # Create Elegant Data Visualisations Using the Grammar of Graphics
  library(gridExtra) # Miscellaneous Functions for "Grid" Graphics
  library(ggpubr) # 'ggplot2' Based Publication Ready Plots
  data$d <- data$T - data$CK
  Fig_Tsum_data <- data.frame(A = c("环境温度(°C)", "增温温度(°C)", "温差(°C)", "t值", "P值"), 
                              B = c(paste(round(mean(data$CK),2), "±", 
                                          round(sd(data$CK)/sqrt(length(data$CK)),2)),
                                    paste(round(mean(data$T),2), "±", 
                                          round(sd(data$T)/sqrt(length(data$T)),2)),
                                    paste(round(mean(data$d),2), "±", 
                                          round(sd(data$d)/sqrt(length(data$d)),2)),
                                    round(t.test(data$CK, data$T, pair=TRUE)[["statistic"]],3),
                                    round(t.test(data$CK, data$T, pair=TRUE)[["p.value"]],3)))
  if (Fig_Tsum_data[5,2] < 0.001) {
    Fig_Tsum_data[5,2] <- "<0.001"
  }
  g_data <- data.frame(date = rep(data$date,2),
                       temperature = c(data$CK, data$T),
                       group = c(rep("环境对照", length(data$date)), 
                                 rep("增温处理", length(data$date))))
  Fig_T <- ggplot(g_data)+
    geom_line(aes(x = date, y = temperature, color = group), size = 1)+
    scale_x_date(#date_breaks = "1 week",
                 date_labels = "%m/%d")+
    scale_y_continuous(limits = c(min(g_data$temperature), 
                                  1.25*max(g_data$temperature) - 0.25*min(g_data$temperature)))+
    labs(x = "", y = "温度(°C)")+
    scale_color_manual(values = color_yb2)+
    annotation_custom(grob = tableGrob(t(Fig_Tsum_data), rows = NULL, cols = NULL),
                      ymin = 1.1*max(g_data$temperature) - 0.1*min(g_data$temperature),
                      ymax = 1.2*max(g_data$temperature) - 0.2*min(g_data$temperature))+
    theme_bw()+
    theme(legend.position="top",
          legend.title = element_blank(),
          legend.text = element_text(size = 12, face = "bold"),
          axis.text.x = element_text(size = 12, face = "bold"),
          axis.title.x = element_blank(),
          axis.title.y = element_text(size = 14),
          axis.text.y = element_text(size = 12),
          text = element_text(family = "myFont"))
  return(Fig_T)
} #温度数据绘图（含差异分析，结果进图）



# Section 4 第一章 综合分析--------------------------------------------------------
# 读取数据
Data1 <- read.xlsx("1-data.xlsx",sheetName = "data", encoding = "UTF-8")  #读取待分析数据
M_ECFs <- data.frame(m = c("DC", "DN", "NC", "TC", "TD", "TN", "PN"),
                     s1 = c("D", "D", "N", "T", "T", "T", "P"),
                     s2 = c("C", "N", "C", "C", "D", "N", "N"))  #两ECFs组合关系

# 统计分析
## 效应值计算：E-Effect size，P-Performance，R-RCI
### 全球变化因子
E_P_ECF <- bootstrap.yb(Data1, "ReTG", "ReTC", group_c = "Target_Type", c("ECFs", "NUM_ECFs_Record"), 1000, 100)
E_C_ECF <- bootstrap.yb(Data1, "TGRCI", "TCRCI", group_c = "Target_Type", c("ECFs", "NUM_ECFs_Record"),  1000, 100)

## 计算差异值Difference size：仅ECF-M与ECF-S之间
Index <- unique(rbindlist(E_P_ECF)[,c(1,3)])
Index$Num <- c(1:nrow(Index))
Index_M <- Index[nchar(Index$ECFs) == 2,]
for (i in c(1:nrow(Index_M))) {
  Index_M$ECF1[i] <- substring(Index_M[i,1], 1, 1)
  Index_M$ECF2[i] <- substring(Index_M[i,1], 2, 2)
}
Index_M$Num1[Index_M$Target_Type == "Invasive"] <- as.numeric(replace.yb(Index_M$ECF1[Index_M$Target_Type == "Invasive"], c("C", "D", "N", "P", "T"), c(1, 3, 9, 13, 17)))
Index_M$Num2[Index_M$Target_Type == "Invasive"] <- as.numeric(replace.yb(Index_M$ECF2[Index_M$Target_Type == "Invasive"], c("C", "D", "N", "P", "T"), c(1, 3, 9, 13, 17)))
Index_M$Num1[Index_M$Target_Type == "Native"] <- as.numeric(replace.yb(Index_M$ECF1[Index_M$Target_Type == "Native"], c("C", "D", "N", "P", "T"), c(2, 4, 10, 14, 18)))
Index_M$Num2[Index_M$Target_Type == "Native"] <- as.numeric(replace.yb(Index_M$ECF2[Index_M$Target_Type == "Native"], c("C", "D", "N", "P", "T"), c(2, 4, 10, 14, 18)))
D_P_ECF <- list()
for (i in c(1:nrow(Index_M))) {
  D_P_ECF[[i]] <- cbind(E_P_ECF[[Index_M$Num[i]]][,c(1:3)],
                        data.frame(D = E_P_ECF[[Index_M$Num[i]]]$E - 
                                     (E_P_ECF[[Index_M$Num1[i]]]$E + E_P_ECF[[Index_M$Num2[i]]]$E)))
}
D_C_ECF <- list()
for (i in c(1:nrow(Index_M))) {
  D_C_ECF[[i]] <- cbind(E_C_ECF[[Index_M$Num[i]]][,c(1:3)],
                        data.frame(D = E_C_ECF[[Index_M$Num[i]]]$E - 
                                     (E_C_ECF[[Index_M$Num1[i]]]$E + E_C_ECF[[Index_M$Num2[i]]]$E)))
}

## 交互判定
I_P_ECF <- inter.yb(E_P_ECF, 
                 c(5, 6, 7, 8, 11, 12, 15, 16, 19, 20, 21, 22, 23, 24), 
                 c(3, 4, 3, 4, 9, 10, 13, 14, 17, 18, 17, 18, 17, 18), 
                 c(1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 3, 4, 9, 10), "E")
I_C_ECF <- inter.yb(E_C_ECF, 
                 c(5, 6, 7, 8, 11, 12, 15, 16, 19, 20, 21, 22, 23, 24), 
                 c(3, 4, 3, 4, 9, 10, 13, 14, 17, 18, 17, 18, 17, 18), 
                 c(1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 3, 4, 9, 10), "E")

## 合并效应值、差异值
SE_P_ECF <- mq.yb(rbindlist(E_P_ECF), "E", c("Target_Type", "ECFs", "NUM_ECFs_Record"))
SE_C_ECF <- mq.yb(rbindlist(E_C_ECF), "E", c("Target_Type", "ECFs", "NUM_ECFs_Record"))
SD_P_ECF <- mq.yb(rbindlist(D_P_ECF), "D", c("Target_Type", "ECFs"))
SD_C_ECF <- mq.yb(rbindlist(D_C_ECF), "D", c("Target_Type", "ECFs"))
SI_P_ECF <- as.data.frame(table(rbindlist(I_P_ECF)$ECFs, rbindlist(I_P_ECF)$Target_Type, rbindlist(I_P_ECF)$inter))
names(SI_P_ECF) <- c("ECFs", "Target_Type", "inter", "Freq")
SI_C_ECF <- as.data.frame(table(rbindlist(I_C_ECF)$ECFs, rbindlist(I_C_ECF)$Target_Type, rbindlist(I_C_ECF)$inter))
names(SI_C_ECF) <- c("ECFs", "Target_Type", "inter", "Freq")
SE_P_F <- mq.yb(rbindlist(E_P_ECF), "E", c("Target_Type", "NUM_ECFs_Record"))
SE_C_F <- mq.yb(rbindlist(E_C_ECF), "E", c("Target_Type", "NUM_ECFs_Record"))
SD_P_F <- mq.yb(rbindlist(D_P_ECF), "D", c("Target_Type", "NUM_ECFs_Record"))
SD_C_F <- mq.yb(rbindlist(D_C_ECF), "D", c("Target_Type", "NUM_ECFs_Record"))

## 获取样本量
N_P_ECF <- ddply(Data1, c("Target_Type", "ECFs"), summarise, N = length(ReTG))
N_C_ECF <- ddply(Data1, c("Target_Type", "ECFs"), summarise, N = length(ReTG))
N_P_F <- ddply(Data1, c("Target_Type", "NUM_ECFs_Record"), summarise, N = length(ReTG))
N_C_F <- ddply(Data1, c("Target_Type", "NUM_ECFs_Record"), summarise, N = length(ReTG))
N_P_N <- ddply(Data1, c("Target_Type", "NFIX_Target"), summarise, N = length(ReTG))
N_C_N <- ddply(Data1, c("Target_Type", "NFIX_Target"), summarise, N = length(ReTG))
N_P_LC <- ddply(Data1, c("Target_Type", "Lifecycle_Target"), summarise, N = length(ReTG))
N_C_LC <- ddply(Data1, c("Target_Type", "Lifecycle_Target"), summarise, N = length(ReTG))
N_P_FG <- ddply(Data1, c("Target_Type", "FUN_Group_Target"), summarise, N = length(ReTG))
N_C_FG <- ddply(Data1, c("Target_Type", "FUN_Group_Target"), summarise, N = length(ReTG))
N_P_I <- ddply(Data1, c("Target_Type", "Type_IND"), summarise, N = length(ReTG))
N_C_I <- ddply(Data1, c("Target_Type", "Type_IND"), summarise, N = length(ReTG))

## 数据整理
P_F <- merge(N_P_F, SE_P_F, by = c("Target_Type", "NUM_ECFs_Record"), all = TRUE) 
C_F <- merge(N_C_F, SE_C_F, by = c("Target_Type", "NUM_ECFs_Record"), all = TRUE) 
P_F$ECFs[P_F$NUM_ECFs_Record == 1] <- "ECF(S)"
P_F$ECFs[P_F$NUM_ECFs_Record == 2] <- "ECF(M)"
C_F$ECFs[C_F$NUM_ECFs_Record == 1] <- "ECF(S)"
C_F$ECFs[C_F$NUM_ECFs_Record == 2] <- "ECF(M)"
P_ECF <- merge(N_P_ECF, SE_P_ECF, by = c("Target_Type", "ECFs"), all = TRUE) 
C_ECF <- merge(N_C_ECF, SE_C_ECF, by = c("Target_Type", "ECFs"), all = TRUE) 
P_ECF <- rbind(merge(N_P_ECF, SE_P_ECF, by = c("Target_Type", "ECFs"), all = TRUE), P_F)
C_ECF <- rbind(merge(N_C_ECF, SE_C_ECF, by = c("Target_Type", "ECFs"), all = TRUE), C_F)
names(SD_P_F)[2] <- "ECFs"
SD_P_F$ECFs[SD_P_F$ECFs == 2] <- "ECF(M)"
SD_P_ECF <- rbind(SD_P_ECF, SD_P_F)
names(SD_C_F)[2] <- "ECFs"
SD_C_F$ECFs[SD_C_F$ECFs == 2] <- "ECF(M)"
SD_C_ECF <- rbind(SD_C_ECF, SD_C_F)


# 绘图
## 图1-3单因子结果
D1.3 <- rbind(
  cbind(P_ECF[P_ECF$NUM_ECFs_Record == 1,], 
        data.frame(PC = rep("(A) Performance", length(P_ECF[P_ECF$NUM_ECFs_Record == 1,1])))),
  cbind(C_ECF[C_ECF$NUM_ECFs_Record == 1,], 
        data.frame(PC = rep("(B) Competitiveness", length(C_ECF[C_ECF$NUM_ECFs_Record == 1,1]))))
)
D1.3$PC <- factor(D1.3$PC, levels = c("(A) Performance", "(B) Competitiveness"))
D1.3$sig[D1.3$p < 0.025 | D1.3$p > 0.975] <- "*"


F1.3 <- ggplot(data = D1.3)+
  geom_point(aes(x = mean, y = ECFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = ECFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, linewidth = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.7, color = "#727272")+
  geom_text(aes(x = up, y = ECFs, label = N, group = Target_Type), color = "#000000",
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(aes(x = down, y = ECFs, label = sig, group = Target_Type), size = 8, color = "#000000",
            position=position_dodge(width = -0.5), hjust = 1.5, vjust = 0.7, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = color_yb2)+
  scale_fill_manual(values = color_yb2)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("ECF(S)","C","N","D","P","T"), 
                                   levels=c("ECF(S)","C","N","D","P","T")))+
  labs(x = "Effect size")+
  facet_wrap(vars(PC),scales= "free_x")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(size=14,),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        strip.text.x = element_text(size=14),
        text = element_text(family = "myFont"))

## 图1.4两因子结果
D1.4 <- rbind(
  cbind(P_ECF[P_ECF$NUM_ECFs_Record == 2,], 
        data.frame(PC = rep("(A) Performance", length(P_ECF[P_ECF$NUM_ECFs_Record == 2,1])))),
  cbind(C_ECF[C_ECF$NUM_ECFs_Record == 2,], 
        data.frame(PC = rep("(B) Competitiveness", length(C_ECF[C_ECF$NUM_ECFs_Record == 2,1]))))
)
D1.4$PC <- factor(D1.4$PC, levels = c("(A) Performance", "(B) Competitiveness"))
D1.4$sig[D1.4$p < 0.025 | D1.4$p > 0.975] <- "*"

F1.4 <- ggplot(data = D1.4)+
  geom_point(aes(x = mean, y = ECFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = ECFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, linewidth = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.7, color = "#727272")+
  geom_text(aes(x = up, y = ECFs, label = N, group = Target_Type), color = "#000000",
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(aes(x = down, y = ECFs, label = sig, group = Target_Type), color = "#000000", size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, vjust = 0.7, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = color_yb2)+
  scale_fill_manual(values = color_yb2)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("ECF(M)","NC","DC","DN","PN","TC","TN","TD"), 
                                   levels=c("ECF(M)","NC","DC","DN","PN","TC","TN","TD")))+
  labs(x = "Effect size")+
  facet_wrap(vars(PC),scales= "free_x")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(size=14,),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        strip.text.x = element_text(size=14),
        text = element_text(family = "myFont"))

## 图1.5交互结果
D1.5 <- rbind(
  cbind(SD_P_ECF, 
        data.frame(PC = rep("(A) Performance", length(P_ECF[P_ECF$NUM_ECFs_Record == 2,1])))),
  cbind(SD_C_ECF, 
        data.frame(PC = rep("(B) Competitiveness", length(C_ECF[C_ECF$NUM_ECFs_Record == 2,1]))))
)
D1.5$PC <- factor(D1.5$PC, levels = c("(A) Performance", "(B) Competitiveness"))
D1.5$sig[D1.5$p < 0.025 | D1.5$p > 0.975] <- "*"
D1.5$N <- D1.4$N

F1.5 <- ggplot(data = D1.5)+
  geom_point(aes(x = mean, y = ECFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = ECFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, linewidth = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.7, color = "#727272")+
  geom_text(aes(x = up, y = ECFs, label = N, group = Target_Type), color = "#000000",
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(aes(x = down, y = ECFs, label = sig, group = Target_Type), color = "#000000", size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, vjust = 0.7, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = color_yb2)+
  scale_fill_manual(values = color_yb2)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("ECF(M)","NC","DC","DN","PN","TC","TN","TD"), 
                                   levels=c("ECF(M)","NC","DC","DN","PN","TC","TN","TD")))+
  labs(x = "Different size")+
  facet_wrap(vars(PC),scales= "free_x")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(size=14,),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        strip.text.x = element_text(size=14),
        text = element_text(family = "myFont"))

# 结果输出
## 表
filename1 <- paste("Results-1 ", Sys.Date(), ".xlsx", sep = "") #第一章结果excel文件名
if (file.exists(filename1)) {
  file.remove(filename1)
}
write.xlsx(P_ECF, filename1, sheetName = "P_ECF")
write.xlsx(C_ECF, filename1, sheetName = "C_ECF", append = TRUE)
write.xlsx(SD_P_ECF, filename1, sheetName = "SD_P_ECF", append = TRUE)
write.xlsx(SD_C_ECF, filename1, sheetName = "SD_C_ECF", append = TRUE)
write.xlsx(SI_P_ECF, filename1, sheetName = "SI_P_ECF", append = TRUE)
write.xlsx(SI_C_ECF, filename1, sheetName = "SI_C_ECF", append = TRUE)

## 图
ggsave(file = "F1.3.png",plot = F1.3, width = 8, height = 6, bg = "white")
ggsave(file = "F1.4.png",plot = F1.4, width = 8, height = 7.5, bg = "white")
ggsave(file = "F1.5.png",plot = F1.5, width = 8, height = 7.5, bg = "white")

# Section 5 第三章 变暖与光周期--------------------------------------------------------

# 数据读取
grow_leaf36 <- read.xlsx("3-data.xlsx",sheetName = "Sheet1", encoding = "UTF-8")  #生长和叶片数据
grow_leaf36$T <- as.factor(grow_leaf36$T)
grow_leaf36$L <- as.factor(grow_leaf36$L)
grow_leaf36$U <- as.factor(grow_leaf36$U)
grow_leaf36$C <- as.factor(grow_leaf36$C)
grow_leaf36$S <- as.factor(grow_leaf36$S)

grow_leaf3 <- grow_leaf36[grow_leaf36$U == 0,]

cnp36 <- read.xlsx("3-data.xlsx",sheetName = "CNP", encoding = "UTF-8")   #元素数据
cnp36$T <- as.factor(cnp36$T)
cnp36$L <- as.factor(cnp36$L)
cnp36$U <- as.factor(cnp36$U)
cnp36$C <- as.factor(cnp36$C)
cnp36$S <- as.factor(cnp36$S)

cnp3 <- cnp36[cnp36$U == 0,]

temperature3 <- read.xlsx("3-Temperature.xlsx", sheetName = "Sheet1")

#计算的数据
grow_leaf3$TB <- grow_leaf3$AB + grow_leaf3$UB  #总生物量
grow_leaf3$RS <- grow_leaf3$UB / grow_leaf3$AB  #根冠比
grow_leaf3$SLA <- (grow_leaf3$LA / 100) / (grow_leaf3$LM) #比叶面积
cnp3$CN <- cnp3$OC / cnp3$TN * 1000
cnp3$CP <- cnp3$OC / cnp3$TP * 1000
cnp3$NP <- cnp3$TN / cnp3$TP

#分析
Results_grow_leaf3 <- analysis.yb(data = grow_leaf3, 
                                  Col_DV = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                                  Col_IV = c("T", "L"), 
                                  Col_group = "C", Col_species = "S", species = c("I", "N"))
Results_cnp3 <- analysis.yb(data = cnp3, 
                            Col_DV = c("CN", "CP", "NP"), 
                            Col_IV = c("T", "L"), 
                            Col_group = "C", Col_species = "S", species = c("I", "N"))


#绘图
Fig3 <- c(fig.yb(data = Results_grow_leaf3[[1]][[1]][[3]], 
                 Col_IV = c("T", "L"), 
                 x1 = c("IV_0_0", "IV_0_1", "IV_0_2", "IV_1_0", "IV_1_1", "IV_1_2"), 
                 x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                 Col_DV = "DV", 
                 y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"),  
                 y2 = c("株高(cm)", "根长(cm)", "基径(mm)", "总生物量(g)", "根冠比", 
                        "叶片数", expression("叶面积 ( " * mm ^ 2 * ")"), "叶干重(g)", 
                        expression("比叶面积 (" * cm ^ 2 * "/g)"), "Fv/Fo", "Fv/Fm", "叶绿素相对含量(SPAD)", "叶氮(mg/g)"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F),
          fig.yb(data = Results_cnp3[[1]][[1]][[3]], 
                 Col_IV = c("T", "L"), 
                 x1 = c("IV_0_0", "IV_0_1", "IV_0_2", "IV_1_0", "IV_1_1", "IV_1_2"), 
                 x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                 Col_DV = "DV", 
                 y1 = c("CN", "CP", "NP"), 
                 y2 = c("C:N", "C:P", "N:P"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F))

Fig_I3 <- c(fig.yb(data = Results_grow_leaf3[[2]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_1", "IV_0_2", "IV_1_0", "IV_1_1", "IV_1_2"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp3[[2]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_1", "IV_0_2", "IV_1_0", "IV_1_1", "IV_1_2"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))

Fig_N3 <- c(fig.yb(data = Results_grow_leaf3[[3]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_1", "IV_0_2", "IV_1_0", "IV_1_1", "IV_1_2"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp3[[3]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_1", "IV_0_2", "IV_1_0", "IV_1_1", "IV_1_2"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))

fig3.1 <- temperature.yb(temperature3) 

fig3.2 <- ggarrange(Fig_I3[[1]], Fig_I3[[2]], Fig_I3[[3]], Fig_I3[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig3.3 <- ggarrange(Fig_N3[[1]], Fig_N3[[2]], Fig_N3[[3]], Fig_N3[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig3.4 <- ggarrange(Fig3[[1]], Fig3[[2]], Fig3[[3]], Fig3[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

fig3.5 <- ggarrange(Fig_I3[[6]], Fig_I3[[7]], Fig_I3[[8]], Fig_I3[[9]], 
                    Fig_I3[[10]], Fig_I3[[11]], Fig_I3[[12]], Fig_I3[[13]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig3.7 <- ggarrange(Fig_N3[[6]], Fig_N3[[7]], Fig_N3[[8]], Fig_N3[[9]], 
                    Fig_N3[[10]], Fig_N3[[11]], Fig_N3[[12]], Fig_N3[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig3.9 <- ggarrange(Fig3[[6]], Fig3[[7]], Fig3[[8]], Fig3[[9]], 
                    Fig3[[10]], Fig3[[11]], Fig3[[12]],  Fig3[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 4, ncol = 2, label.y = 1.05)

fig3.10 <- ggarrange(Fig_I3[[5]], Fig_I3[[14]], Fig_I3[[15]], Fig_I3[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig3.11 <- ggarrange(Fig_N3[[5]], Fig_N3[[14]], Fig_N3[[15]], Fig_N3[[16]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig3.12 <- ggarrange(Fig3[[5]], Fig3[[14]], Fig3[[15]], Fig3[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

# 结果输出
## 表
filename3 <- paste("Results-3 ", Sys.Date(), ".xlsx", sep = "")
sheetname3 <- data.frame(c("F", "Sum"),
                         c("I-F", "I-Sum"),
                         c("N-F", "N-Sum"))

if (file.exists(filename3)) {
  file.remove(filename3)
}
for (i in 1:3) {
  for (j in 2:3) {
    if (i == 1 & j == 1) {
      write.xlsx(rbind(Results_grow_leaf3[[i]][[1]][[j]], Results_cnp3[[i]][[1]][[j]]),
                 file = filename3, sheetName = as.character(sheetname3[j-1,i]), append = FALSE, row.names = FALSE)
    } else {
      write.xlsx(rbind(Results_grow_leaf3[[i]][[1]][[j]], Results_cnp3[[i]][[1]][[j]]),
                 file = filename3, sheetName = as.character(sheetname3[j-1,i]), 
                 append = TRUE, row.names = FALSE)
    }
  }
}

## 图
fig_set3 <- data.frame(name = c("fig3.1", "fig3.2", "fig3.3", "fig3.4", "fig3.5", 
                                "fig3.7", "fig3.9", "fig3.10", "fig3.11", "fig3.12"),
                      width = c(12,15,15,15,15,15,15,15,15,15),
                      height = c(4,3,3,6,6,6,12,3,3,6))
for (i in c(1:10)) {
  ggsave(file = paste(fig_set3$name[i], ".png", sep = ""), 
         plot = eval(parse(text = fig_set3$name[i])), 
         width = fig_set3$width[i], height = fig_set3$height[i],  bg = "white")
}


# Section 6 第四章 变暖与光强--------------------------------------------------------

grow_leaf45 <- read.xlsx("5-data.xlsx",sheetName = "Sheet1", encoding = "UTF-8")  #生长和叶片数据
grow_leaf45$T <- as.factor(grow_leaf45$T)
grow_leaf45$L <- as.factor(grow_leaf45$L)
grow_leaf45$C <- as.factor(grow_leaf45$C)
grow_leaf45$S <- as.factor(grow_leaf45$S)

grow_leaf4 <- grow_leaf45[grow_leaf45$L %in% c(0,1,2),]

cnp45 <- read.xlsx("5-data.xlsx",sheetName = "CNP", encoding = "UTF-8")   #元素数据
cnp45$T <- as.factor(cnp45$T)
cnp45$L <- as.factor(cnp45$L)
cnp45$C <- as.factor(cnp45$C)
cnp45$S <- as.factor(cnp45$S)

cnp4 <- cnp45[cnp45$L %in% c(0,1,2),]

temperature4 <- read.xlsx("5-Temperature.xlsx", sheetName = "Sheet1")

#计算的数据
grow_leaf4$TB <- grow_leaf4$AB + grow_leaf4$UB  #总生物量
grow_leaf4$RS <- grow_leaf4$UB / grow_leaf4$AB  #根冠比
grow_leaf4$SLA <- (grow_leaf4$LA / 100) / (grow_leaf4$LM) #比叶面积
cnp4$CN <- cnp4$OC / cnp4$TN * 1000
cnp4$CP <- cnp4$OC / cnp4$TP * 1000
cnp4$NP <- cnp4$TN / cnp4$TP

#分析
Results_grow_leaf4 <- analysis.yb(data = grow_leaf4, 
                                  Col_DV = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                                  Col_IV = c("T", "L"), 
                                  Col_group = "C", Col_species = "S", species = c("I", "N"))
Results_cnp4 <- analysis.yb(data = cnp4, 
                            Col_DV = c("CN", "CP", "NP"), 
                            Col_IV = c("T", "L"), 
                            Col_group = "C", Col_species = "S", species = c("I", "N"))

#结果输出图片
Fig4 <- c(fig.yb(data = Results_grow_leaf4[[1]][[1]][[3]], 
                 Col_IV = c("T", "L"), 
                 x1 = c("IV_0_1", "IV_0_0", "IV_0_2", "IV_1_1", "IV_1_0", "IV_1_2"), 
                 x2 = c("T0LL","T0LN",  "T0LH", "T1LL", "T1LN", "T1LH"),
                 Col_DV = "DV", 
                 y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"),  
                 y2 = c("株高(cm)", "根长(cm)", "基径(mm)", "总生物量(g)", "根冠比", 
                        "叶片数", expression("叶面积 (" * mm ^ 2 * ")"), "叶干重(g)", 
                        expression("比叶面积 (" * cm ^ 2 * "/g)"), "Fv/Fo", "Fv/Fm", "叶绿素相对含量(SPAD)", "叶氮(mg/g)"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F),
          fig.yb(data = Results_cnp4[[1]][[1]][[3]], 
                 Col_IV = c("T", "L"), 
                 x1 = c("IV_0_1", "IV_0_0", "IV_0_2", "IV_1_1", "IV_1_0", "IV_1_2"), 
                 x2 = c("T0LL","T0LN",  "T0LH", "T1LL", "T1LN", "T1LH"), 
                 Col_DV = "DV", 
                 y1 = c("CN", "CP", "NP"), 
                 y2 = c("C:N", "C:P", "N:P"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F))

Fig_I4 <- c(fig.yb(data = Results_grow_leaf4[[2]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_1", "IV_0_0", "IV_0_2", "IV_1_1", "IV_1_0", "IV_1_2"), 
                   x2 = c("T0LL","T0LN",  "T0LH", "T1LL", "T1LN", "T1LH"),
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp4[[2]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_1", "IV_0_0", "IV_0_2", "IV_1_1", "IV_1_0", "IV_1_2"), 
                   x2 = c("T0LL","T0LN",  "T0LH", "T1LL", "T1LN", "T1LH"),
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))

Fig_N4 <- c(fig.yb(data = Results_grow_leaf4[[3]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_1", "IV_0_0", "IV_0_2", "IV_1_1", "IV_1_0", "IV_1_2"), 
                   x2 = c("T0LL","T0LN",  "T0LH", "T1LL", "T1LN", "T1LH"),
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp4[[3]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_1", "IV_0_0", "IV_0_2", "IV_1_1", "IV_1_0", "IV_1_2"), 
                   x2 = c("T0LL","T0LN",  "T0LH", "T1LL", "T1LN", "T1LH"),
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))

fig4.1 <- temperature.yb(temperature4) 

fig4.2 <- ggarrange(Fig_I4[[1]], Fig_I4[[2]], Fig_I4[[3]], Fig_I4[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig4.3 <- ggarrange(Fig_N4[[1]], Fig_N4[[2]], Fig_N4[[3]], Fig_N4[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig4.4 <- ggarrange(Fig4[[1]], Fig4[[2]], Fig4[[3]], Fig4[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

fig4.5 <- ggarrange(Fig_I4[[6]], Fig_I4[[7]], Fig_I4[[8]], Fig_I4[[9]], Fig_I4[[10]], Fig_I4[[11]], Fig_I4[[12]], Fig_I4[[13]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig4.6 <- ggarrange(Fig_N4[[6]], Fig_N4[[7]], Fig_N4[[8]], Fig_N4[[9]], Fig_N4[[10]], Fig_N4[[11]], Fig_N4[[12]], Fig_N4[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig4.7 <- ggarrange(Fig4[[6]], Fig4[[7]], Fig4[[8]], Fig4[[9]], Fig4[[10]], Fig4[[11]], Fig4[[12]],  Fig4[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 4, ncol = 2, label.y = 1.05)

fig4.8 <- ggarrange(Fig_I4[[5]], Fig_I4[[14]], Fig_I4[[15]], Fig_I4[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig4.9 <- ggarrange(Fig_N4[[5]], Fig_N4[[14]], Fig_N4[[15]], Fig_N4[[16]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig4.10 <- ggarrange(Fig4[[5]], Fig4[[14]], Fig4[[15]], Fig4[[16]], 
                     align = "hv", labels = LETTERS, font.label = list(size = 16),
                     common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

# 结果输出
## 表
filename4 <- paste("Results-4 ", Sys.Date(), ".xlsx", sep = "")
if (file.exists(filename4)) {
  file.remove(filename4)
}
for (i in 1:3) {
  for (j in 2:3) {
    if (i == 1 & j == 1) {
      write.xlsx(rbind(Results_grow_leaf4[[i]][[1]][[j]], Results_cnp4[[i]][[1]][[j]]),
                 file = filename4, sheetName = as.character(sheetname3[j-1,i]), 
                 append = FALSE, row.names = FALSE)
    } else {
      write.xlsx(rbind(Results_grow_leaf4[[i]][[1]][[j]], Results_cnp4[[i]][[1]][[j]]),
                 file = filename4, sheetName = as.character(sheetname3[j-1,i]), 
                 append = TRUE, row.names = FALSE)
    }
  }
}

## 图
fig_set4 <- data.frame(name = c("fig4.1", "fig4.2", "fig4.3", "fig4.4", "fig4.5", 
                                "fig4.6", "fig4.7", "fig4.8", "fig4.9", "fig4.10"),
                       width = c(12,15,15,15,15,15,15,15,15,15),
                       height = c(4,3,3,6,6,6,12,3,3,6))
for (i in c(1:10)) {
  ggsave(file = paste(fig_set4$name[i], ".png", sep = ""), 
         plot = eval(parse(text = fig_set4$name[i])), 
         width = fig_set4$width[i], height = fig_set4$height[i],  bg = "white")
}


# Section 7 第五章 变暖与光质--------------------------------------------------------

grow_leaf5 <- grow_leaf45[grow_leaf45$L %in% c(0,3,4),]
cnp5 <- cnp45[cnp45$L %in% c(0,3,4),]

light <- read.xlsx("6-light.xlsx",sheetName = "Sheet2", encoding = "UTF-8")  #生长和叶片数据
light2 <- data.frame()
for (i in c(1:30)) {
  T <- rep(light[1,i+1],401)
  L <- rep(light[2,i+1],401)
  re <- rep(light[3,i+1],401)
  wave <- light$NA.[c(-1,-2,-3)]
  int <- light[c(-1,-2,-3),i+1]
  light2 <- rbind(light2, data.frame(T, L, re, wave, int))
  rm(T, L, re, wave, int)
}
light2$wave <- as.numeric(light2$wave)
light2$int <- as.numeric(light2$int)

#计算的数据
grow_leaf5$TB <- grow_leaf5$AB + grow_leaf5$UB  #总生物量
grow_leaf5$RS <- grow_leaf5$UB / grow_leaf5$AB  #根冠比
grow_leaf5$SLA <- (grow_leaf5$LA / 100) / (grow_leaf5$LM) #比叶面积
cnp5$CN <- cnp5$OC / cnp5$TN * 1000
cnp5$CP <- cnp5$OC / cnp5$TP * 1000
cnp5$NP <- cnp5$TN / cnp5$TP

#分析
Results_grow_leaf5 <- analysis.yb(data = grow_leaf5, 
                                  Col_DV = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                                  Col_IV = c("T", "L"), 
                                  Col_group = "C", Col_species = "S", species = c("I", "N"))
Results_cnp5 <- analysis.yb(data = cnp5, 
                            Col_DV = c("CN", "CP", "NP"), 
                            Col_IV = c("T", "L"), 
                            Col_group = "C", Col_species = "S", species = c("I", "N"))


# 绘图
Fig5 <- c(fig.yb(data = Results_grow_leaf5[[1]][[1]][[3]], 
                 Col_IV = c("T", "L"), 
                 x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                 x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                 Col_DV = "DV", 
                 y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                 y2 = c("株高(cm)", "根长(cm)", "基径(mm)", "总生物量(g)", "根冠比", 
                        "叶片数", expression("叶面积 (" * mm ^ 2 * ")"), "叶干重(g)", 
                        expression("比叶面积 (" * cm ^ 2 * "/g)"), "Fv/Fo", "Fv/Fm", "叶绿素相对含量(SPAD)", "叶氮(mg/g)"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F),
          fig.yb(data = Results_cnp5[[1]][[1]][[3]], 
                 Col_IV = c("T", "L"), 
                 x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                 x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                 Col_DV = "DV", 
                 y1 = c("CN", "CP", "NP"), 
                 y2 = c("C:N", "C:P", "N:P"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F))

Fig_I5 <- c(fig.yb(data = Results_grow_leaf5[[2]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp5[[2]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))

Fig_N5 <- c(fig.yb(data = Results_grow_leaf5[[3]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp5[[3]][[1]][[3]], 
                   Col_IV = c("T", "L"), 
                   x1 = c("IV_0_0", "IV_0_3", "IV_0_4", "IV_1_0", "IV_1_3", "IV_1_4"), 
                   x2 = c("T0L0", "T0L1", "T0L2", "T1L0", "T1L1", "T1L2"), 
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))

light_fig <- ddply(light2, .(L, wave), summarize, int = mean(int))
light_fig <- light_fig[light_fig$L %in% c(0,3,4),]
light_fig$L[light_fig$L == 0] <- "L0 (CK)"
light_fig$L[light_fig$L == 3] <- "L1 (R)"
light_fig$L[light_fig$L == 4] <- "L2 (FR)"

fig5.1 <- ggplot(light_fig)+
  geom_col(aes(x = wave, y = int, fill = wave, colour = wave), show.legend = FALSE)+
  scale_fill_gradientn(values = seq(0,1,0.05),
                       colours = c('#610061','#8300b5','#6a00ff','#0000ff','#007bff','#00d5ff','#00ff92',
                                   '#36ff00','#81ff00','#c3ff00','#ffff00','#ffbe00','#ff7700','#ff2100',
                                   '#ff0000','#ff0000','#ff0000','#db0000','#b50000','#8d0000','#610000'))+
  scale_color_gradientn(values = seq(0,1,0.05),
                        colours = c('#610061','#8300b5','#6a00ff','#0000ff','#007bff','#00d5ff','#00ff92',
                                    '#36ff00','#81ff00','#c3ff00','#ffff00','#ffbe00','#ff7700','#ff2100',
                                    '#ff0000','#ff0000','#ff0000','#db0000','#b50000','#8d0000','#610000'))+
  facet_grid(L ~ .)+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0), limits =c(0,390))+
  labs(x = "波长(nm)", y = expression("光照强度(μW/  " * cm ^ 2 * "/nm)"))+
  theme_bw()+
  theme(axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        strip.text = element_text(size=12, face="bold"))

fig5.2 <- ggarrange(Fig_I5[[1]], Fig_I5[[2]], Fig_I5[[3]], Fig_I5[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig5.3 <- ggarrange(Fig_N5[[1]], Fig_N5[[2]], Fig_N5[[3]], Fig_N5[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig5.4 <- ggarrange(Fig5[[1]], Fig5[[2]], Fig5[[3]], Fig5[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

fig5.5 <- ggarrange(Fig_I5[[6]], Fig_I5[[7]], Fig_I5[[8]], 
                    Fig_I5[[9]], Fig_I5[[10]], Fig_I5[[11]], Fig_I5[[12]], Fig_I5[[13]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig5.6 <- ggarrange(Fig_N5[[6]], Fig_N5[[7]], Fig_N5[[8]], 
                    Fig_N5[[9]], Fig_N5[[10]], Fig_N5[[11]], Fig_N5[[12]], Fig_N5[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig5.7 <- ggarrange(Fig5[[6]], Fig5[[7]], Fig5[[8]], Fig5[[9]], 
                    Fig5[[10]], Fig5[[11]], Fig5[[12]],  Fig5[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 4, ncol = 2, label.y = 1.05)

fig5.8 <- ggarrange(Fig_I5[[5]], Fig_I5[[14]], Fig_I5[[15]], Fig_I5[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig5.9 <- ggarrange(Fig_N5[[5]], Fig_N5[[14]], Fig_N5[[15]], Fig_N5[[16]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig5.10 <- ggarrange(Fig5[[5]], Fig5[[14]], Fig5[[15]], Fig5[[16]], 
                     align = "hv", labels = LETTERS, font.label = list(size = 16),
                     common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

# 结果输出
## 表
filename5 <- paste("Results-5 ", Sys.Date(), ".xlsx", sep = "")
if (file.exists(filename5)) {
  file.remove(filename5)
}
for (i in 1:3) {
  for (j in 2:3) {
    if (i == 1 & j == 1) {
      write.xlsx(rbind(Results_grow_leaf5[[i]][[1]][[j]], Results_cnp5[[i]][[1]][[j]]),
                 file = filename5, sheetName = as.character(sheetname3[j-1,i]), 
                 append = FALSE, row.names = FALSE)
    } else {
      write.xlsx(rbind(Results_grow_leaf5[[i]][[1]][[j]], Results_cnp5[[i]][[1]][[j]]),
                 file = filename5, sheetName = as.character(sheetname3[j-1,i]), 
                 append = TRUE, row.names = FALSE)
    }
  }
}

## 图
fig_set5 <- data.frame(name = c("fig5.1", "fig5.2", "fig5.3", "fig5.4", "fig5.5", 
                                "fig5.6", "fig5.7", "fig5.8", "fig5.9", "fig5.10"),
                       width = c(10,15,15,15,15,15,15,15,15,15),
                       height = c(7,3,3,6,6,6,12,3,3,6))
for (i in c(1:10)) {
  ggsave(file = paste(fig_set5$name[i], ".png", sep = ""), 
         plot = eval(parse(text = fig_set5$name[i])), 
         width = fig_set5$width[i], height = fig_set5$height[i],  bg = "white")
}







# Section 8 第六章 变暖与UV-A--------------------------------------------------------

# 数据读取
grow_leaf6 <- grow_leaf36[grow_leaf36$L == 0,]
cnp6 <- cnp36[cnp36$L == 0,]

# 计算的数据
grow_leaf6$TB <- grow_leaf6$AB + grow_leaf6$UB  #总生物量
grow_leaf6$RS <- grow_leaf6$UB / grow_leaf6$AB  #根冠比
grow_leaf6$SLA <- (grow_leaf6$LA / 100) / (grow_leaf6$LM) #比叶面积
cnp6$CN <- cnp6$OC / cnp6$TN * 1000
cnp6$CP <- cnp6$OC / cnp6$TP * 1000
cnp6$NP <- cnp6$TN / cnp6$TP

# 统计分析
Results_grow_leaf6 <- analysis.yb(data = grow_leaf6, 
                                  Col_DV = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                                  Col_IV = c("T", "U"), 
                                  Col_group = "C", Col_species = "S", species = c("I", "N"))
Results_cnp6 <- analysis.yb(data = cnp6, 
                            Col_DV = c("CN", "CP", "NP"), 
                            Col_IV = c("T", "U"), 
                            Col_group = "C", Col_species = "S", species = c("I", "N"))


#结果输出图片
Fig6 <- c(fig.yb(data = Results_grow_leaf6[[1]][[1]][[3]], 
                 Col_IV = c("T", "U"), 
                 x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                 x2 = c("CK", "T", "U", "T×U"), 
                 Col_DV = "DV", 
                 y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                 y2 = c("株高(cm)", "根长(cm)", "基径(mm)", "总生物量(g)", "根冠比", 
                        "叶片数", expression("叶面积 (" * mm ^ 2 * ")"), "叶干重(g)", 
                        expression("比叶面积 (" * cm ^ 2 * "/g)"), "Fv/Fo", "Fv/Fm", "叶绿素相对含量(SPAD)", "叶氮(mg/g)"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F),
          fig.yb(data = Results_cnp6[[1]][[1]][[3]], 
                 Col_IV = c("T", "U"), 
                 x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                 x2 = c("CK", "T", "U", "T×U"), 
                 Col_DV = "DV", 
                 y1 = c("CN", "CP", "NP"), 
                 y2 = c("C:N", "C:P", "N:P"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_I","Mix_N", "Mono_I", "Mono_N"), 
                 g2 = c("混合种植-鬼针草","混合种植-金盏银盘", "单一种植-鬼针草", "单一种植-金盏银盘"), 
                 pointrange = F))

Fig_I6 <- c(fig.yb(data = Results_grow_leaf6[[2]][[1]][[3]], 
                   Col_IV = c("T", "U"), 
                   x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                   x2 = c("CK", "T", "U", "T×U"), 
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp6[[2]][[1]][[3]], 
                   Col_IV = c("T", "U"), 
                   x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                   x2 = c("CK", "T", "U", "T×U"), 
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))

Fig_N6 <- c(fig.yb(data = Results_grow_leaf6[[3]][[1]][[3]], 
                   Col_IV = c("T", "U"), 
                   x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                   x2 = c("CK", "T", "U", "T×U"), 
                   Col_DV = "DV", 
                   y1 = c("Ht", "RL", "BS", "TB", "RS", "LN", "LA", "LM", "SLA", "Fv.Fo", "Fv.Fm", "Chl", "ChlN"), 
                   y2 = c("株高", "根长", "基径", "总生物量", "根冠比", 
                          "叶片数", "叶面积", "叶干重", 
                          "比叶面积", "Fv/Fo", "Fv/Fm", "叶绿素相对含量", "叶氮"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T),
            fig.yb(data = Results_cnp6[[3]][[1]][[3]], 
                   Col_IV = c("T", "U"), 
                   x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                   x2 = c("CK", "T", "U", "T×U"), 
                   Col_DV = "DV", 
                   y1 = c("CN", "CP", "NP"), 
                   y2 = c("C:N", "C:P", "N:P"), 
                   Col_group = c("C"), 
                   g1 = c("Mix", "Mono"), 
                   g2 = c("混合种植", "单一种植"), 
                   pointrange = T))


fig6.1 <- ggarrange(Fig_I6[[1]], Fig_I6[[2]], Fig_I6[[3]], Fig_I6[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.2 <- ggarrange(Fig_N6[[1]], Fig_N6[[2]], Fig_N6[[3]], Fig_N6[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.3 <- ggarrange(Fig6[[1]], Fig6[[2]], Fig6[[3]], Fig6[[4]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

fig6.4 <- ggarrange(Fig_I6[[6]], Fig_I6[[7]], Fig_I6[[8]], Fig_I6[[9]], 
                    Fig_I6[[10]], Fig_I6[[11]], Fig_I6[[12]], Fig_I6[[13]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig6.5 <- ggarrange(Fig_N6[[6]], Fig_N6[[7]], Fig_N6[[8]], Fig_N6[[9]], 
                    Fig_N6[[10]], Fig_N6[[11]], Fig_N6[[12]], Fig_N6[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig6.6 <- ggarrange(Fig6[[6]], Fig6[[7]], Fig6[[8]], Fig6[[9]], 
                    Fig6[[10]], Fig6[[11]], Fig6[[12]],  Fig6[[13]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 4, ncol = 2, label.y = 1.05)

fig6.7 <- ggarrange(Fig_I6[[5]], Fig_I6[[14]], Fig_I6[[15]], Fig_I6[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.8 <- ggarrange(Fig_N6[[5]], Fig_N6[[14]], Fig_N6[[15]], Fig_N6[[16]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig6.9 <- ggarrange(Fig6[[5]], Fig6[[14]], Fig6[[15]], Fig6[[16]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

# 结果输出
## 表
filename6 <- paste("Results-6 ", Sys.Date(), ".xlsx", sep = "")

if (file.exists(filename6)) {
  file.remove(filename6)
}
for (i in 1:3) {
  for (j in 2:3) {
    if (i == 1 & j == 1) {
      write.xlsx(rbind(Results_grow_leaf6[[i]][[1]][[j]], Results_cnp6[[i]][[1]][[j]]),
                 file = filename6, sheetName = as.character(sheetname3[j-1,i]), append = FALSE, row.names = FALSE)
    } else {
      write.xlsx(rbind(Results_grow_leaf6[[i]][[1]][[j]], Results_cnp6[[i]][[1]][[j]]),
                 file = filename6, sheetName = as.character(sheetname3[j-1,i]), 
                 append = TRUE, row.names = FALSE)
    }
  }
}

## 图
fig_set6 <- data.frame(name = c("fig6.1", "fig6.2", "fig6.3", "fig6.4", 
                                "fig6.5", "fig6.6", "fig6.7", "fig6.8", "fig6.9"),
                      width = c(12,12,12,12,12,12,12,12,12),
                      height = c(3,3,6,6,6,12,3,3,6))
for (i in c(1:9)) {
  ggsave(file = paste(fig_set6$name[i], ".png", sep = ""), 
         plot = eval(parse(text = fig_set6$name[i])), 
         width = fig_set6$width[i], height = fig_set6$height[i],  bg = "white")
}


# Section 9 第七章 变暖与UV-B--------------------------------------------------------

# 数据读取
grow7 <- read.xlsx("2-data.xlsx",sheetName = "grow", encoding = "UTF-8")  #生长数据
grow7$T <- as.factor(grow7$T)
grow7$U <- as.factor(grow7$U)
grow7$C <- as.factor(grow7$C)
grow7$S <- as.factor(grow7$S)

leaf7 <- read.xlsx("2-data.xlsx",sheetName = "leaf", encoding = "UTF-8")  #叶片数据
leaf7$T <- as.factor(leaf7$T)
leaf7$U <- as.factor(leaf7$U)
leaf7$C <- as.factor(leaf7$C)
leaf7$S <- as.factor(leaf7$S)

ps7 <- read.xlsx("2-data.xlsx",sheetName = "ps", encoding = "UTF-8")     #光合数据
ps7$T <- as.factor(ps7$T)
ps7$U <- as.factor(ps7$U)
ps7$C <- as.factor(ps7$C)
ps7$S <- as.factor(ps7$S)

cnp7 <- read.xlsx("2-data.xlsx",sheetName = "cnp", encoding = "UTF-8")   #元素数据
cnp7$T <- as.factor(cnp7$T)
cnp7$U <- as.factor(cnp7$U)
cnp7$C <- as.factor(cnp7$C)
cnp7$S <- as.factor(cnp7$S)

temperature7 <- read.xlsx("2-Temperature.xlsx", sheetName = "Sheet1")

# 计算的数据
grow7$TB <- grow7$AB + grow7$UB #总生物量
grow7$RS <- grow7$UB / grow7$AB #根冠比
cnp7$CN <- cnp7$OC / cnp7$TN * 1000
cnp7$CP <- cnp7$OC / cnp7$TP * 1000
cnp7$NP <- cnp7$TN / cnp7$TP

# 统计分析
Results_grow7 <- analysis.yb(data = grow7, 
                             Col_DV = c("H", "RL", "BS", "TB", "RS"), 
                             Col_IV = c("T", "U"), 
                             Col_group = "C", Col_species = "S", species = c("SC", "AA"))
Results_leaf7 <- analysis.yb(data = leaf7, 
                             Col_DV = c("leaf_a", "leaf_m"), 
                             Col_IV = c("T", "U"), 
                             Col_group = "C", Col_species = "S", species = c("SC", "AA"))
Results_ps7 <- analysis.yb(data = ps7, 
                           Col_DV = c("Fv.Fm", "A", "Ci", "E", "gs", "WUE", "chl_N", "chl"), 
                           Col_IV = c("T", "U"), 
                           Col_group = "C", Col_species = "S", species = c("SC", "AA"))
Results_cnp7 <- analysis.yb(data = cnp7, 
                            Col_DV = c("CN", "CP", "NP"), 
                            Col_IV = c("T", "U"), 
                            Col_group = "C", Col_species = "S", species = c("SC", "AA"))


# 绘图

Fig7 <- c(fig.yb(data = Results_grow7[[1]][[1]][[3]], 
                 Col_IV = c("T", "U"), 
                 x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                 x2 = c("CK", "T", "U", "T×U"), 
                 Col_DV = "DV", 
                 y1 = c("H", "RL", "BS", "TB", "RS"), 
                 y2 = c("株高(cm)", "根长(cm)", "基径(mm)", "总生物量(mg)", "根冠比"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_SC","Mix_AA", "Mono_SC", "Mono_AA"), 
                 g2 = c("混合种植-加拿大一枝黄花","混合种植-艾草", "单一种植-加拿大一枝黄花", "单一种植-艾草"), 
                 pointrange = F),
          fig.yb(data = Results_leaf7[[1]][[1]][[3]], 
                 Col_IV = c("T", "U"), 
                 x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                 x2 = c("CK", "T", "U", "T×U"), 
                 Col_DV = "DV", 
                 y1 = c("leaf_a", "leaf_m"), 
                 y2 = c(expression("叶面积 ( " * mm ^ 2 * ")"), "叶鲜重(mg)"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_SC","Mix_AA", "Mono_SC", "Mono_AA"), 
                 g2 = c("混合种植-加拿大一枝黄花","混合种植-艾草", "单一种植-加拿大一枝黄花", "单一种植-艾草"), 
                 pointrange = F),
          fig.yb(data = Results_ps7[[1]][[1]][[3]], 
                 Col_IV = c("T", "U"), 
                 x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                 x2 = c("CK", "T", "U", "T×U"), 
                 Col_DV = "DV", 
                 y1 = c("Fv.Fm", "A", "Ci", "E", "gs", "WUE", "chl_N", "chl"), 
                 y2 = c("Fv/Fm", 
                        expression("净光合速率 ( " * μmol ~ m ^ -2 ~ s ^ -1 * ")"),
                        expression("胞间 " * CO[2] *" 浓度(ppm)"), 
                        expression("蒸腾速率 ( " * mmol ~ m ^ -2 ~ s ^ -1 * ")"),
                        expression("气孔导度 ( " * mmol ~ m ^ -2 ~ s ^ -1 * ")"),
                        "水分利用效率(%)", "叶氮(mg/g)", "叶绿素相对含量(SPAD)"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_SC","Mix_AA", "Mono_SC", "Mono_AA"), 
                 g2 = c("混合种植-加拿大一枝黄花","混合种植-艾草", "单一种植-加拿大一枝黄花", "单一种植-艾草"), 
                 pointrange = F),
          fig.yb(data = Results_cnp7[[1]][[1]][[3]], 
                 Col_IV = c("T", "U"), 
                 x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                 x2 = c("CK", "T", "U", "T×U"), 
                 Col_DV = "DV", 
                 y1 = c("CN", "CP", "NP"), 
                 y2 = c("C:N", "C:P", "N:P"), 
                 Col_group = c("C", "S"), 
                 g1 = c("Mix_SC","Mix_AA", "Mono_SC", "Mono_AA"), 
                 g2 = c("混合种植-加拿大一枝黄花","混合种植-艾草", "单一种植-加拿大一枝黄花", "单一种植-艾草"), 
                 pointrange = F))

Fig_SC7 <- c(fig.yb(data = Results_grow7[[2]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("H", "RL", "BS", "TB", "RS"), 
                    y2 = c("株高", "根长", "基径", "总生物量", "根冠比"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T),
             fig.yb(data = Results_leaf7[[2]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("leaf_a", "leaf_m"), 
                    y2 = c("叶面积", "叶鲜重"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T),
             fig.yb(data = Results_ps7[[2]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("Fv.Fm", "A", "Ci", "E", "gs", "WUE", "chl_N", "chl"), 
                    y2 = c("Fv/Fm", "净光合速率", expression("胞间 " * CO[2] *"浓度"), 
                           "蒸腾速率", "气孔导度", "水分利用效率", "叶氮", "叶绿素相对含量"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T),
             fig.yb(data = Results_cnp7[[2]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("CN", "CP", "NP"), 
                    y2 = c("C:N", "C:P", "N:P"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T))

Fig_AA7 <- c(fig.yb(data = Results_grow7[[3]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("H", "RL", "BS", "TB", "RS"), 
                    y2 = c("株高", "根长", "基径", "总生物量", "根冠比"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T),
             fig.yb(data = Results_leaf7[[3]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("leaf_a", "leaf_m"), 
                    y2 = c("叶面积", "叶鲜重"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T),
             fig.yb(data = Results_ps7[[3]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("Fv.Fm", "A", "Ci", "E", "gs", "WUE", "chl_N", "chl"), 
                    y2 = c("Fv/Fm", "净光合速率", expression("胞间 " * CO[2] *"浓度"), 
                           "蒸腾速率", "气孔导度", "水分利用效率", "叶氮", "叶绿素相对含量"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T),
             fig.yb(data = Results_cnp7[[3]][[1]][[3]], 
                    Col_IV = c("T", "U"), 
                    x1 = c("IV_0_0", "IV_1_0", "IV_0_1", "IV_1_1"), 
                    x2 = c("CK", "T", "U", "T×U"), 
                    Col_DV = "DV", 
                    y1 = c("CN", "CP", "NP"), 
                    y2 = c("C:N", "C:P", "N:P"), 
                    Col_group = c("C"), 
                    g1 = c("Mix", "Mono"), 
                    g2 = c("混合种植", "单一种植"), 
                    pointrange = T))

fig7.1 <- temperature.yb(temperature7) 

fig7.2 <- ggarrange(Fig_SC7[[1]], Fig_SC7[[2]], Fig_SC7[[3]], 
                    Fig_SC7[[4]], Fig_SC7[[6]], Fig_SC7[[7]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 3, label.y = 1.05)
fig7.3 <- ggarrange(Fig_AA7[[1]], Fig_AA7[[2]], Fig_AA7[[3]], 
                    Fig_AA7[[4]], Fig_AA7[[6]], Fig_AA7[[7]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = TRUE, legend = "top", nrow = 2, ncol = 3, label.y = 1.05)
fig7.4 <- ggarrange(Fig7[[1]], Fig7[[2]], Fig7[[3]], 
                    Fig7[[4]], Fig7[[6]], Fig7[[7]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = TRUE, nrow = 3, ncol = 2, label.y = 1.05)

fig7.5 <- ggarrange(Fig_SC7[[8]], Fig_SC7[[9]], Fig_SC7[[10]], 
                    Fig_SC7[[11]], Fig_SC7[[12]], Fig_SC7[[13]], 
                    Fig_SC7[[14]], Fig_SC7[[15]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig7.6 <- ggarrange(Fig_AA7[[8]], Fig_AA7[[9]], Fig_AA7[[10]], 
                    Fig_AA7[[11]], Fig_AA7[[12]], Fig_AA7[[13]], 
                    Fig_AA7[[14]], Fig_AA7[[15]],  
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 2, ncol = 4, label.y = 1.05)
fig7.7 <- ggarrange(Fig7[[8]], Fig7[[9]], Fig7[[10]], 
                    Fig7[[11]], Fig7[[12]], Fig7[[13]], 
                    Fig7[[14]], Fig7[[15]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 4, ncol = 2, label.y = 1.05)

fig7.8 <- ggarrange(Fig_SC7[[5]], Fig_SC7[[16]], Fig_SC7[[17]], Fig_SC7[[18]], 
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig7.9 <- ggarrange(Fig_AA7[[5]], Fig_AA7[[16]], Fig_AA7[[17]], Fig_AA7[[18]],
                    align = "hv", labels = LETTERS, font.label = list(size = 16),
                    common.legend = T, legend = "top", nrow = 1, ncol = 4, label.y = 1.05)
fig7.10 <- ggarrange(Fig7[[5]], Fig7[[16]], Fig7[[17]], Fig7[[18]],
                     align = "hv", labels = LETTERS, font.label = list(size = 16),
                     common.legend = T, legend = "top", nrow = 2, ncol = 2, label.y = 1.05)

# 结果输出
## 表
filename7 <- paste("Results-7 ", Sys.Date(), ".xlsx", sep = "")
if (file.exists(filename7)) {
  file.remove(filename7)
}
sheetname7 <- data.frame(c("F", "Sum"),
                         c("SC-F", "SC-Sum"),
                         c("AA-F", "AA-Sum"))
for (i in 1:3) {
  for (j in 2:3) {
    if (i == 1 & j == 1) {
      write.xlsx(rbind(Results_grow7[[i]][[1]][[j]], Results_leaf7[[i]][[1]][[j]], 
                       Results_ps7[[i]][[1]][[j]], Results_cnp7[[i]][[1]][[j]]),
                 file = filename7, sheetName = as.character(sheetname7[j-1,i]), 
                 append = FALSE, row.names = FALSE)
    } else {
      write.xlsx(rbind(Results_grow7[[i]][[1]][[j]], Results_leaf7[[i]][[1]][[j]], 
                       Results_ps7[[i]][[1]][[j]], Results_cnp7[[i]][[1]][[j]]),
                 file = filename7, sheetName = as.character(sheetname7[j-1,i]), 
                 append = TRUE, row.names = FALSE)
    }
  }
}

## 图
fig_set7 <- data.frame(name = c("fig7.1", "fig7.2", "fig7.3", "fig7.4", "fig7.5", 
                                "fig7.6", "fig7.7", "fig7.8", "fig7.9", "fig7.10"),
                       width = c(12,9,9,12,12,12,12,12,12,12),
                       height = c(4,6,6,9,6,6,12,3,3,6))
for (i in c(1:10)) {
  ggsave(file = paste(fig_set7$name[i], ".png", sep = ""), 
         plot = eval(parse(text = fig_set7$name[i])), 
         width = fig_set7$width[i], height = fig_set7$height[i],  bg = "white")
}

